import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
from scipy.stats import uniform
from numpy import log,exp,pi

# Initialize random number generator
np.random.seed(123)

# True parameter values
mu_x, mu_y = 0, 0
sigma_x, sigma_y = 1, 1
rho_true = 0.5

# Generate data
N=1000
data=np.random.multivariate_normal([mu_x,mu_y],[[sigma_x, rho_true],[rho_true, sigma_y]], N)
x=data[:,0]
y=data[:,1]

E=10000 # Number of iteration
BURN_IN=50 # First number of samples to be droped 

# Initialize the chain. 
rho=0

# Store the samples
chain_rho=np.array([0.]*(E-BURN_IN))

accepted_number=0.
for e in range(E):
    # Draw a value from the proposal distribution, Uniform(rho-0.07,rho+0.07); Equation 7
    rho_candidate=uniform.rvs(rho-0.1,2*0.1)

    # Compute the acceptance probability, Equation 8 and Equation 6. 
    # We will do both equations in log domain here to avoid underflow.
    accept=-rho_candidate**2- N*log((1.-rho_candidate**2)**(1./2)) - sum(1./(2.*(1.-rho_candidate**2))*(x**2-2.*rho_candidate*x*y+y**2))
    accept=accept-(-rho**2- N*log((1.-rho**2)**(1./2)) - sum(1./(2.*(1.-rho**2))*(x**2-2.*rho*x*y+y**2)))
    accept=min([0,accept])
    accept=exp(accept)

    # Accept rho_candidate with probability accept.
    if uniform.rvs(0,1)<accept:
        rho=rho_candidate
        accepted_number=accepted_number+1
    else:
        rho=rho

    # store
    if e>=BURN_IN:
        chain_rho[e-BURN_IN]=rho

print("Acceptance ratio is "+str(accepted_number/(E)))
print("Mean rho is "+str(chain_rho.mean()))
print("Std for rho is "+str(chain_rho.std()))
print("Compare with np.cov function: "+str(np.cov(data.T)))

Acceptance ratio is 0.3322
Mean rho is 0.5234692091037678
Std for rho is 0.020250036657334457
Compare with np.cov function: [[0.98010969 0.50207627]
 [0.50207627 0.96368557]]

fig = plt.figure(figsize=(10,4))
ax1 = fig.add_subplot(311)
ax1.scatter(x,y,s=20,c='b',marker = 'o')
ax2 = fig.add_subplot(312)
ax2.plot(chain_rho,'b')
ax2.set_ylabel('$rho$')
ax3 = fig.add_subplot(313)
ax3.hist(chain_rho,50)
ax3.set_xlabel('$rho$')
plt.show()